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1 Introduction 



1.1 QCD and Jet Fractions 

QCD is now well accepted as the theory of strong interactions. It has a running coupling 
constant as — o^s{Q 2 )/2tt, which is only low enough for perturbation theory to be valid 
at high energies: Q 2 is the relevant momentum scale. In the LEP and SLC e + e~ colliders, 
the centre-of-mass energy Q > 90 GeV is high enough for perturbation theory to be valid 
in the early stages of the strong interaction [|IJ: it can therefore be used to calculate the 
production of quarks and gluons (partons). These energies also allow the determination of 
the QCD colour factors Cp and Ca 0, [|- The partons then combine at lower momentum 
scales where perturbation theory is not valid into a large number of hadrons in the process 
known as hadronisation. It is these hadrons that are actually measured in the detectors. 

The distribution of experimentally measured hadrons closely follows the theoretically cal- 
culated partons in its flow of momentum and quantum numbers |4], |5|. The uncertainty 
principle relates the distance over which an interaction takes place to the momentum trans- 
fer involved. Interactions with a momentum scale greater than about 1 GeV are in the 
perturbative regime. In the non-perturbative regime there are large distances between 
partons, and the momentum transfer at this stage is therefore small. This phenomenon, 
known as Local Parton-Hadron Duality, allows hadronisation to be modelled by Monte 
Carlo event generators such as JETSET || [7| and HERWIG ^ [], |HJ . The former uses a 
string model in which the gluon field energy is converted into hadrons, whilst the latter uses 
a low-mass cluster model. These two event generators yield similar results, for example in 



the determination of as [O, [L2L [L3|, [L4| , giving some confidence in their validity. 



One of the common observables used in the study of hadronic final states is the number of 
jets as defined by a given jet algorithm. The nth jet fraction R n is then defined as 

Rn = — (1) 

O total 

where a n is the cross section for n jets and a tota i is the total cross section. Jet fractions may 
be used to calculate as by comparing parton-level theoretical calculations with measured 
hadron-level rates, provided the same jet algorithm is used in the two stages of the analysis. 

Jet algorithms use a measure to determine the closeness of pairs of tracks. The JADE 
algorithm [[T^, [Tj|] uses the invariant mass 



j = (Pi+Pj) 2 = 2^ (?] 
Uij Q 2 Q 9 1 ' 



where Pi and pj are the four-momenta of the two particles (assumed massless) and Q 2 is 
the centre-of-mass energy squared of the complete event as defined above. Let the partons 
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have energies E{, Ej and relative angle 9 in the centre-of-mass frame. The measure is then 



j _ 2E l E J {l- C os6) 

Vij - q 2 ■ V 6 ) 

At each stage, the pair of tracks with the smallest measure are combined if this measure 
is smaller than a fixed cut-off, ie if < y cut . The algorithm stops once all pairs have 
a measure greater than y cut , at which point the remaining tracks are termed jets. In 
experimental analyses the tracks are hadrons in the detector; in theoretical calculations 
they are partons. The validity of comparing these two situations is at least suggested 
by Local Parton-Hadron Duality, and might provide information about non-perturbative 
effects. 



The Durham algorithm |17j is identical to the JADE algorithm except that it uses as its 
measure 



D 2min{^ 2 ,fff}(l-cos^ 



A further variant on this theme is the Geneva algorithm JE8|, again identical except for its 
measure: 

G = 8 £^(1 - 008 0) 

Uij 9 (Et + E,f ' 1 ] 

This algorithm has not been widely used in experimental analyses. For each measure there 
are a number of different possible procedures for combining pairs of tracks, the simplest 
of which is to add their four-momenta. The calculations presented in this paper are of an 
accuracy not to be affected by the procedure used. 

The mass of a current quark (~ few MeV) is much less than the characteristic QCD 
scale (Aqcd — 200 MeV). For the above algorithms the effect of the masses is swamped 
by non-perturbative effects fl2~0 |, and the quarks may be considered massless. 



At small values of y cut the jet fractions are dominated by soft gluons, which have a useful 
property known as angular ordering. Consider the emission of such a soft gluon from 
either of a pair of partons separated by an angle 9. The emitting partons may be treated 
as classical currents: the emission from each of these currents interferes, but is coherent 



due to the colour structure. The resulting cross section pi] , written in terms of the 
energies and angles of the partons, may be split into two parts, each containing a collinear 
singularity. When averaged over the azimuthal direction the contribution from each of 
these parts is non-zero only within a cone of angle 9: the emission may be separated into 
two parts, each of which describes emission from a single parton, but restricted to the cone 
described above. Emission at larger angles adds up to be equivalent to emission from the 
parent parton ||22|| . This argument holds to next-to-leading log order, and was discovered 



independently by Mueller [23fl and Ermolaev and Fadin . Coherence effects of this type 



have been observed at a variety of energies [25, 26, 27 . 
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Figure 1: Jet rates as measured by the OPAL collaboration. Data taken from 128 



1.2 The Leading Logarithmic Approximation 

Hadronisation corrections set in at y cut ~ m2 H a dronslQ 2 ~ (1 GeV) 2 /Q 2 : at LEP energies 
this gives a lower limit on y cut of ~ 0.0001. At values of y cut < 0.1, such that 

L = In ( — J > 2.5, (6) 

V Vcut J 

the divergent regions of phase space due to soft, collinear gluons dominate the integrand. 
Perturbation theory then gives a series in two variables, in increasing powers of ag and 
decreasing powers of L. For jet fractions the leading term in L has two powers of L for 
every power of as, one from each of the soft and collinear divergences. At each order of as 
the series in L can be truncated to give a valid approximation, provided that asL 2 < 1. 
I shall work to leading log order, in which only the highest power of L is kept. At LEP 



energies as — 0.117 ± 0.005 ||29|| , giving a lower limit on y cut of ~ 0.05. At higher energies 



this limit will be lowered, widening the valid range of y cut . 

At lower values of y cu t at LEP energies it is necessary to resum the series to all orders 
in ag- This is possible if the coefficients follow, or are at least close to, the pattern of 
a known function; for example the leading log coefficients for the two jet fraction in the 
Durham algorithm follow the pattern of the coefficients of the exponential function p0| . 
This phenomenon is known as exponentiation pH] , and depends on the factorisation of 
the phase space of the relevant observable: the QCD matrix elements are guaranteed to 
factorise for soft, collinear gluons These resummed predictions compare well with 
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data p3[, and may be used to extract the value of as, the main error being theoretical. 



The same properties do not hold for the JADE algorithm, which has more complicated 
phase space boundaries due to extra gluon-gluon correlations: the coefficients follow no 
discernable pattern for the most commonly used recombination prescriptions [34j]. 



Calculations for jet fractions exist at leading order in as and leading log order up to six jet 



order for the JADE algorithm |p5| , and up to four jet order for the Durham algorithm |L7 , 



36]. The JADE calculations only include the abelian terms, that is only the diagrams with 
no triple gluon vertices. The existing results are here extended to include the non-abelian 
terms. The existence of the triple gluon vertex has been confirmed at LEP [[|. Some of 
the coefficients have also been investigated numerically [[STf] . 

Jet clustering has also been used as a first stage in experimental analyses to reduce the 



theoretical errors in other observables pq| . The calculation of jet fractions can also be used 
as a test of the consistency of QCD between different experimental situations, for example 
between e + e~ annihilation and deep inelastic scattering. 



1.3 Soft and Collinear Divergences 

In the perturbative calculation of parton branching there are two types of divergence that 
occur even aside from renormalisation. In the case of branching from an outgoing parton 
(time-like branching), and after integrating over the azimuthal direction, the cross section 
can be expressed as [ |39| 

da n+1 _ sin 0d9 

— = a s V{z)dz- (7) 

aa n 1 — cos 9 

where 9 is the angle between the partons. For q — > qg branching z is the energy fraction of 
the emitted gluon, whilst for g — > gg branching it is the energy fraction of either gluon. The 
formulation in terms of parton energies and angles is used in order to take into account both 
types of divergence. V(z) is the splitting function which gives the probability of parton 
branching after averaging and summing over initial and final spin states respectively. In 
the soft (small z) limit the two splitting functions relevant to this analysis [^TJ are 



V?(s)[? - qg] = C F 1 + (1 "" )2 * ^ + 0(1) 



and 



V§{z)[g^gg] = 2C A 



z(l-z) 



l-z 



1 - z 



^ + 0(1). (9) 



These splitting functions diverge at z = 0. The evolution equation (^) is also divergent in 
the collinear (small 9) limit: 

dcr n+i Q _ T/ , , , d0 dzd9 2 

— ~ 2a s V{z)dz— ~ 2a s C — , 10 

da„ 9 z 9 l 
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where C = Cf or Ca according to the type of branching. This divergence is a consequence 
of the partons' masslessness, but the expression would be logarithmically enhanced even 
with small parton masses. The splitting function for gluon to quark- ant iquark branching 
is 

VZ{z)\g-+qq\ = \N f [# + (l-zf\, (11) 

where Nf is the number of active quark flavours at the relevant energy. This function 
contains no soft divergence, and therefore this type of branching does not contribute at 
leading log order. Similarly the four gluon vertex is suppressed by a factor of as- 

For an observable to be physical it must not be affected by any of the divergences in the 
theory. It must therefore take the same value if any parton is replaced by an arbitrary 
number of partons such that the total momentum and colour flow of the group remains the 
same. Such an observable is called collinear-safe. For the same reason it must also remain 
unchanged if a parton emits an arbitrary number of gluons of zero energy, in which case 
it is called infrared- safe. This suggests that a parton is significant for its momentum and 
colour flow, rather than as a directly measurable particle. 



2 Analytic Calculation of Jet Fractions 

For the leading log order result all the gluons are soft and collinear. The soft limit enables 
the energy fraction of each parton to be approximated as e« ~ Ei/Q. Note that the 
original quark and antiquark both have e = 1 to lowest order in e. For two partons 
branching off the same initial parton the collinear limit (9 small) allows the approximation 
2(1 — cos 9) ~ 9 2 + 0{9 4 ) to be made. The measures (|3[) and (|) then simplify to 

VZ = (12) 
yg = mm{e 2 ,e 2 }9 2 . (13) 

For two gluons originating one from the quark and one from the antiquark, the collinear 
limit is 9 ~ 7r and the measures become 

Vtj = (14) 



yB = min{ej, ej}. (15) 



P. _ ...... 

These approximations are sufficient for the leading log order calculation. In the small 9 
limit the perpendicular momentum of the less energetic particle with respect to the more 
energetic particle is k± = Qe sin 9 ~ Qe9 where Qe is the energy of the less energetic par- 
ticle as defined above. The Durham measure is thus seen to be equal to the perpendicular 
momentum squared divided by Q 2 . 

The phase space for n jets fllCf ) is 

dei de n d9\ d9 2 n 

ir-Tir-ir (16) 
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The substitutions 




XiL, 




y%L, In 



Ucut 



(17) 



transform this to 



dx\...dx n dy\...dy n . (18) 

The formal integration range for the original variables 9f is taken as — > 1, remembering 
that only small values of e and 9 contribute at leading log order. This translates into Xi, yi 
space as —>■ oo. In the actual calculation it is simplest to put in these limits as step 
functions, ie Q(xi)...Q(x n )Q(yi)...Q(y n ), where 



x > 0, 
x < 0. 



(19) 



I will work for the moment with the JADE algorithm. Consider two gluons radiated 
successively by a quark with energy fractions q and €j, and at angles 9i and 9j respectively 
to the quark. In order to generate the maximum number of logarithms L the angles must 
be strongly ordered, ie 9j <C 9,i [21]. The smaller of the angles (9j) may therefore be 
ignored, and the angle between the two partons taken as 0j. The condition that these two 
gluons are resolved as separate jets is y^ > y cut , which appears in the integrand as 



6(e i e j 6 l J 2 - y cut ); 



the substitutions ([L7]) transform this to 



6(1 - Xi - xj - y^. 



(20) 



(21) 



This divides the phase space into two halves separated by a plane. Taken together with 
the appropriate conditions for other pairs of partons the overall integration region is a 
'polygon' in a space with two dimensions for every emitted gluon. 

For the specific diagram 
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the conditions that no two partons will be combined by the JADE algorithm are 



angular ordering 


0(02 - 9 3 ) 


quark-gluon 




gluon-gluon 


©(eie 2 - Vcut) 


0(eie 3 - y cut ) 


©(e 2 e 3 ^3 - Vcut) 



This gives as the integration 

/oo 
0(2/3 - 2/2)0(1 - xx - yi)0(l - an - x 2 )0(l - a:i - ^3)6(1 - x 2 - x 3 - y 2 ) 
-00 

Q(x 1 )<d(x 2 )Q(x 3 )Q(y 1 )0(y2)Q(y3)dx 1 dx2dx 3 dy 1 dy 2 dy 3 . (22) 

This calculation may be done using the C program outlined in appendix (]S|): the answer 
is 53/360, which yields a contribution of 

h 3 C 2 F C A —a 3 s L 6 = — (a s L 2 ) 3 (23) 
2 F 360 5 135 V ; V ' 

to the final answer. The factor of 1/2 is statistical, and results from gluons 2 and 3 being 
identical. In this context a pair of gluons are identical if the swapping of their coordinates 
leaves the integrand invariant: such gluons must either contribute a statistical factor of 
1/2 or have their energies ordered. To reach physically meaningful results it is necessary 
to add up the results from each of the diagrams contributing at the relevant order. At 
leading log order for n jets this consists of all the tree level diagrams involving the original 
quark- ant iquark pair and n — 2 emitted gluons; the creation of further quark-antiquark 
pairs is suppressed by one factor of L as there is no associated soft divergence. 

It is of course possible to integrate over the variables in any order. The accuracy of the 
program had been tested for a number of integrands: it gives the same answer for several 
different integration orders. The results presented in this paper have been found using at 
least two different integration orders, always with consistent results. All the results have 
also been checked numerically, using the VEGAS Monte Carlo routine in Fortran 77. Here 
the integration was performed directly in e, 9 space. Again the results are consistent within 
the errors of the Monte Carlo. 



3 Moment Generating Function Method 

The moment generating function for a probability distribution p(n) is defined as 

00 

<j>{u) = 5>>(n). (24) 

n=0 
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The probabilities themselves are recovered using the formula 



p(n) = -r 
n! 



du n 



(#«)) 



(25) 



u=0 



For the case of a gg pair produced in e + e~ annihilation, the moment generating function 
for the Durham algorithm may be calculated to next-to-leading order using the coherent 
branching formalism [ 40| : it is <p 2 (Q,u) = u 



2 e 2F (Q,«) ) where n 



d i f Q 

Fi(Q,u) = — / r q (Q,q)[(f> g (q,u) 
au JQo 



r,(Q,?) 



Here 



r„(Q,<z) 



4CVa;s 
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In 



<f> g (q,u) 



Q 



- 1] dq 

dq, i > 1. 



(26) 



(27) 



is the q ^ qg branching probability, and <f) g (q,u) is the moment generating function for 
a single gluon. At leading log order as may be taken as a constant, evaluated at the 
momentum scale of the interaction. Defining 



Gi(Q,u) 



d i 
du l 



u 



it may easily be calculated that 



Go 
G 1 
G 2 
G 3 
G^ 

Gr 



u 2 ] e 2F ° 



2 ] 

2u + 2u 2 F 1 ] e 2Fo 

2 + 8uF 1 + u 2 (2F 2 + AF 2 )] e 2F ° 

12F X + u {\2F 2 + 2AF 2 ) + u 2 (2F 3 + 12F 1 F 2 + 8F?)] e 2F ° 
24F 2 + A8F 2 + u (I6F3 + 96F!F 2 + 64F\ 3 ) 
+ u 2 (2F 4 + I6F1F3 + 12F 2 + A8F 2 F 2 + 16^)] e 2F ° 
[40F 3 + 240FxF 2 + 160^ + u (20F 4 + I6OF1F3 + 120F; 
+ A80F 2 F 2 + 160F*) + O (u 2 )] e 2F ° 
[6OF4 + 44OF1F3 + 360F| + 1200F 2 F 2 + 320F? + O (u)] e 



?2 



2F 



(29) 



where the and Fj are evaluated at (Q,u). The quark Sudakov form factor is calculated 

as 

3 F (Q,0) = - J§ o r q (Q,q)dq 



A ? (Q) 



(30) 



using (f) g (q, 0) = as calculated in the next section (f41|) . The jet rates themselves are found 
using equation (|2~5"D to be 



R 2 {Q) = A 2 q (Q) 
R 3 (Q) = (2F 1 )A 2 q (Q) 
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Ra(Q) = (F 2 + 2F 2 )A 2 (Q) 

R*(Q) = Qf 3 + 2F 1 f 2 + |f 1 3 )a;(Q) 

R 6 (Q) = (^F 4 + ^F 1 F 3 + lFi + \F*F 2 + |F X 4 ) A^Q), 



where the Fi are evaluated at (Q, 0). 



The gluon moment generating function is [17 



,(Q, u) = [u + u 2 D (Q, u)] e 



C (Q,u) 



where 



Q 



r/(g)- 



dq 



Qo 



r/(g) 



Qo 



du 1 



dq, 



Ei{Q,u) 



0_ 

du 1 



Qo 



{[2T q (Q, q) - T g {Q, q)} [<j> g (q, u)-l]+T f (q)} dq 



[2T q (Q,q)-T 9 (Q,q)} 



Qo 



dq, i > 1, 



and 



Q(Q,u) 



d i 
du 1 



Qo 



{T g (Q,q) [<f> g (q,u) -1]-T f (q)}dq 



Qo 

r ff (Q,<z) 



Here 



and 



r,(Q,?) 



Q 



c?g, z > 1. 
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r/(g) 



2N f a s 
3g 



are g ^ gg and g —>■ qq branching probabilities respectively. Defining 



Ai(g,«) = ^- [0 9 (g,w)] 



it may easily be calculated that 



A (Q,«) 
Ai(Q,u) 

A 3 (Q,u) 



u 



+ u 2 Aj 

1 + u (Ci + 2D ) + u 2 (Dx + AA)] e Co 
2d + 2A> + u (C 2 + C 2 + 4Di + 4AA) 
+ u 2 (,D 2 + 2D 1 Ci + A)C 2 + AjC?)] e C ° 
[3C 2 + 3C 2 + 6D1 + 6D d + u (C 3 + 3CiC 2 + Cf + 6D 2 
+ 12/^ + 6AA + QD C 2 ) + (w 2 )] e c ° 
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Figure 2: Quantity dependencies in the moment generating function method 



MQ, u ) = + 12CiC 2 + + 12D 2 + 2AD l C l + 12 AA + 12AA* 

+ 0(u)]e c °, (39) 

where the Aj, Cj and Dj are evaluated at (Q, «)• The gluon Sudakov form factor is 
calculated as 

A g (Q) = e c ^ = e -S^ g {Q, q ) + v f{q )d<i (40) 

again using <p g (q, 0) = 0. The Ai are evaluated at (Q, 0) to be 

A o (Q,0) = 
A(Q,0) = A fl (Q) 
A 2 (Q,0) = 2(d + £)o)A ff (Q) 



4»(Q,0) = 3(<7 2 + Cf + 2L> 1 + 2AA)A fl (Q) 
A 4 (Q,0) = 4(C 3 + 3C 1 C 2 + C* + 3D 2 + 6D 1 C 1 + 3D C 2 + 3D C*)A g (Q), (41) 

where the Cj and _Dj are also evaluated at (Q,0). Defining a further form factor as 

A f (Q) = e E °( Q ' 0) = e - lQ 2r i(Q'i)- r !>(Q^- T f(i) d i = A g^ ? (42) 

A 5 (Q) 

and defining 

fl*(Q,«) = ^[e^ w - B) ] I (43) 
it may similarly be calculated that 

B (Q,0) = A f (Q) 
B 1 (Q,0) = (E 1 )A f (Q) 

B 2 (Q,0) = (E 2 + Ef)A f (Q), (44) 
where the Ei are also evaluated at (Q, 0). 

Figure (§) shows the dependencies of each of the quantities Ai to E{. The quantities may 
be calculated cyclically in the order in which they are given. The key point is that A n 
only depends (explicitly or implicitly) on A m -i at u — 0, and so the recursive nature 
of equations ( |3"2"D to (^) is avoided. It is interesting to note how the colour factors are 
introduced by the different quantities: for example each factor of Ci(Q,0) introduces a 
colour factor of Ca through its dependency on T g (Q, q) (equation 
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4 Results 



C F and Ca are the QCD colour factors, as = ds/^Ti and L = \n(l/y cut ). and R% 
are the leading log order n jet rates for the JADE and Durham algorithms respectively. 
The JADE results were found by direct integration as described in section (Q), and the 
Durham results using both that and the moment generating function method described 
in section @: the two methods agree at all calculated orders. In addition, as mentioned 
earlier, the results for both algorithms have been checked numerically using the Monte 
Carlo Integrator VEGAS, again with agreement within the errors of the Monte Carlo. 



Rl 



i?4 



R 5 



Ri 



D 



R 



R 



D 



R 



D 



R 



D 



{C F a s L 2 ) 



1 (Ca 
2 + 6 \C F 



(C F a s L 2 Y 



31 5[Ca\ J_(Ca 
45 18 \C F J 36 \C F 



2\3 



(C F a s L 2 ) 



571 4091 
+ 



2520 20160 \C F J 11520 \C F 



+ 



643 (C 



+ 



107 (G 



20160 V C F 



(C F a s L 2 y 



(C F a s L 2 ) 



l_ (Ca 
2 ^ 12 \C F 

6 12 \C F 
1 1 (C A 



(C F a s L 2 Y 
' 90 \C F 



2\3 



(C F a s L 2 ) 



+ 



7 (C 



17 



cr 



(C F a s L 2 )\ 



(45) 



24 24 \C F J 480 VCf/ 10080 \C F/ 
An additional result is the abelian contribution to the 7 jet fraction for the two algorithms: 

1093 



Rj 



R® 



18900 
1 



120 



(C F a s L 2 ) 5 
{C F a s L 2 f. 



(46) 
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Whilst these abelian contributions cannot be used on their own, they provide further 
evidence that the JADE coefficients do not follow an obvious pattern, while the Durham 
abelian coefficients exponentiate. 

One qualitative aspect of the Durham algorithm is that soft, collinear gluons are combined 
only such that at least one of the pair has no other partons closer to it in angle than its 
partner. This suggests that the phase space boundaries given by comparing other pairs of 
gluons should have no effect on the result: this was confirmed for all the results quoted 
above. 

These results are in complete agreement with the results of Brown and Stirling as presented 



in [35j and [|36| for the JADE and Durham algorithms respectively, and also with the results 
of |L7| for the Durham algorithm, there calculated using moment generating functions as 
in section (Q). The results presented here include, in addition to these previous results, the 
non-abelian (Ca) terms for the JADE algorithm, the five and six jet terms for the Durham 
algorithm, and the abelian seven jet terms for both algorithms. The abelian results, the 
integrands for which are generated automatically, are limited by the memory available to 
the computer. The integrands for the non-abelian results are put in by hand, and the 
results are therefore limited by the number of relevant diagrams: the seven jet calculation 
involves 32 such diagrams, compared with only 13 for six jets. 

For the specific case of SU(3), Cp = 4/3 and Ca = 3. The results of equation (^5|) then 
become 

R J 2 ~ 1.00 

R( ~ 2.67 (a s L 2 ) 

R( ~ 3.33 (a s L 2 ) 2 

R( ~ 3.45 (a s L 2 ) 3 

Rl ~ 3.24 (a s L 2 ) 4 

R° ~ 1.00 

Rf ~ 1.33 (a s L 2 ) 

B% ~ 1.22 (a s L 2 ) 2 

R° ~ 0.97 (a s L 2 ) 3 

Rg ~ 0.72 (a s L 2 )\ (47) 

Again there is no clear pattern in the JADE results. Although the numerical coefficients 
do appear to be under control, it is difficult to see how they might be resummed. The root 
of the problem lies in the necessity to compare all gluon pairings in the definition of the 
phase space. Here the Durham algorithm is an improvement, in that fewer pairs of gluons 
need to be compared. 

One further new result is the sequence of coefficients for a cascade generated by a single 
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gluon in the Durham algorithm: 



R? 9 



1 



R? 9 




r° 9 





226800 (° A(XsL2 ) ■ 



(48) 



The coefficients decrease less rapidly than the exponential series owing to the increasing 
number of triple gluon vertices. 

5 Conclusion 

The jet fractions for the JADE and Durham algorithms have been calculated up to six 
jets in the leading order in as and leading log order approximation. The coefficients for 
the JADE algorithm exhibit no easily discernable pattern, suggesting that they may be 
difficult to resum. 

I am most grateful to G. P. Salam and B. R. Webber for many useful discussions. I would 
also like to thank the CERN Theory Division for hospitality while part of this work was 
carried out. 
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A Program Details 



A.l Step Functions and Recursive Branching 

Where it appears as part of an integrand, a step function is defined mathematically as a 
change in the limits of the integration, ie. 

/oo /»oo 
G(x-a)f(x)dx= / f(x)dx. (49) 
oo J a 

A more complicated situation arises when there are two step functions, either of which 
may give the lower limit: 



j — c 



^( X -a)e( X - mx )d^{pW X x (50) 
The above result can be written in a combined form as 

OO /»00 /'OO 

e(x-a)e(x-b)f(x)dx = e(a-b) / f(x)dx + G(b-a) / f(x)dx. (51) 

-oo J a Jb 

In an integration over more than one dimension the constants a and b can become functions 
of the remaining variables of integration. The information already assumed about a and b 
is written in terms of step functions and the integration over the next dimension is therefore 
of the same form. An integration in n dimensions has been branched recursively to two 
n — 1 dimensional integrations. 

The scheme works analogously for upper limits: 



J — c 



Q(a-x)Q(b-x)f(x)dx= { b (52) 

J_oo J\ x )dx b < a. 



If both limits are given by step functions, the step function giving the lower limit must 
switch on before the step function giving the upper limit switches off: 

/oo rb 
G(x-a)e(b-x)f(x)dx = e(b-a) / f(x)dx. (53) 
-oo J a 

The most general case involves an unspecified number of both upper and lower limits. The 
answer contains one branch for every possible pairing of an individual upper limit and an 
individual lower limit. For each branch there are three sources of step functions in addition 
to any step functions contained implicitly in f(x): 

• the lower limit must be greater than other possible lower limits, 

• the upper limit must be smaller than other possible upper limits, 

• the upper limit must be greater than the lower limit. 
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A. 2 Linked Lists 



At each stage of the integration it is not practical to predict the number of terms that 
the polynomial component of the integrand will contain, nor the number of step functions 
needed. Therefore it is essential to have a data structure that allows for an arbitrary 
number of terms at each stage, for example a linked list. The individual terms are stored 



separately, each one with pointers to the next and previous terms. C |41| , f42fl supports 
dynamic memory allocation: functions are included which can insert or remove terms at 
any desired location, memory being taken from or returned to the pool of free memory as 
required. The recursive design of the integration ensures that only one integrand need be 
stored for each dimension at a time, thus minimising the memory use. 



A. 3 Fractions and Euclid's Algorithm 

The initial integrand for each calculation consists numerically purely of integers; real num- 
bers do not occur. This is the case for both the JADE and Durham algorithms. The 
integrand can therefore be expressed purely using rational numbers at all stages of the 
calculation, including the final answer. The handling of fractions is done by a separate 
program module; fractions are cancelled using Euclid's algorithm for the highest common 



factor (HCF) of two integers |43|, f|4| , as described below: 



Let the two numbers whose HCF is to be found be x\ and x 2 , such that x\ > x 2 . The 
procedure is best described schematically: 

Xi = m 1 X2 + x 3 
x 2 = m 2 x 3 + x 4 

x n = m n x n+1 . (54) 

At each stage the first number is divided by the second number, the remainder being given 
by the third number. Any common factor of x\ and x 2 must divide X\ — m,\x 2 = x%, and 
is hence a common factor of x 2 and x 3 . Any common factor of x 2 and £3 must also be a 
factor of xi, and hence the HCF of x\ and x 2 is the same as the HCF of x 2 and £3. This 
argument is repeated for all equations in the above series to show that the original HCF 
required is also the HCF of x n and x n +i : this is clearly x n+ \ itself. 
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